/*
 * case.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_PAM_s1_v2_g67.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_PAM_s1_v2_g67
{
debug deb;
global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

#define Sqrt sqrt
#define Power pow

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double n;
	double t;
	double tPrim;
	double U;
	double V;
	double h;
	double epsilonF;
	double beta;

	unsigned int L;
};

struct variables {
	double lnF;
	//double ln;
	//double lmC;
	double lmF;
	double nF;			// ok
	double mu;			// ok
	double mF;			// ok
	double dF;			// ok

    double Energy;
	double magnetization;			// ok
 };

/*
 *********************************************************
 *********************** ANALOG **************************
 *********************************************************
*/
/*
volatile double   n = 1.60;
volatile double   t = -1.;
volatile double   tPrim = 0.;
volatile double   U = 9.;
volatile double   V = -2.;
volatile double   h = 0.1;
volatile double   epsilonF = -4.;
volatile double   beta = 500;
volatile int      L = 256;

volatile double   lnF = -3.7948168;
volatile double   lmF = 0.3664242;
volatile double   ln = 0;
volatile double   lmC = 0;
volatile double   nF = 0.8749191;
volatile double   mu = -0.5637287;
volatile double   mF = 0.2149817;
volatile double   dF = 0.2032936;

volatile double def_n;
volatile double def_m;
volatile double Energy;

double fq(int s) {
   return (2.*Power(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))),2))/((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double fdqdnF(int s) {
   return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
           (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/
                 Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
             1.*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
             2.*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
         (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double fdqdmF(int s) {
   return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
             Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
           (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
             1.*s*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
             2.*s*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
                Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
         (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double fdqddF(int s) {
   return 0. + (4.*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
         (0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
           (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/
            (2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
       ((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double fE(double epsilonk, int s, int l) {
   return (epsilonk - 2*ln - lnF - 2*mu + epsilonF - (lmC + lmF + 2*h)*s -
           l*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)))/2.;
}

double fdEdlnF(double epsilonk, int s, int l) {
   return (-(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s)) -
           Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)))/
         (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdlmC(double epsilonk, int s, int l) {
   return (s*(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s) -
             Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2))))/
         (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdlmF(double epsilonk, int s, int l) {
   return -(s*(l*(epsilonk + lnF - epsilonF - lmC*s + lmF*s) +
           Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2))))/
      (2.*Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}

double fdEdq(double epsilonk, int s, int l) {
   return -((l*Power(V,2))/Sqrt(4*Power(V,2)*fq(s) + Power(epsilonk + lnF - epsilonF - lmC*s + lmF*s,2)));
}
*/


/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

inline void clear(parameters &p) {
	p.n = p.t = p.tPrim = p.U = p.V = p.h = p.epsilonF = p.beta = 0.; p.L = 0;
}

inline void clear(variables &v) {
	v.lnF = v.lmF = v.nF = v.mu = v.mF = v.dF = 0;
	v.Energy = v.magnetization = 0;
}

inline void addTo(const variables &q, variables &v, double d=1.) {
	v.lnF += q.lnF*d;
	v.lmF += q.lmF*d;
	v.nF += q.nF*d;
	v.mu += q.mu*d;
	v.mF += q.mF*d;
	v.dF += q.dF*d;
	v.Energy += q.Energy*d;
	v.magnetization += q.magnetization*d;
}

inline void copyTo(variables &v, const variables &q, double d=1.) {
	v.lnF = q.lnF*d;
	v.lmF = q.lmF*d;
	v.nF = q.nF*d;
	v.mu = q.mu*d;
	v.mF = q.mF*d;
	v.dF = q.dF*d;
	v.Energy = q.Energy*d;
	v.magnetization = q.magnetization*d;
}

inline void copyTo(parameters &v, const parameters &q) {
	v.n = q.n;
	v.t = q.t;
	v.tPrim = q.tPrim;
	v.U = q.U;
	v.V = q.V;
	v.h = q.h;
	v.epsilonF = q.epsilonF;
	v.beta = q.beta;
	v.L = q.L;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters  p.n  p.t  p.tPrim  p.U  p.V  p.h  p.epsilonF  p.beta  p.L
	variables v; // variables  v.lnF  v.lnF  v.lmC  v.lmF  v.nF  v.mu  v.mC  v.mF  v.d  v.Energy v.magnetization
	double nF, dF, mF, mu, lnF, lmF; // contractions

	// Abbreviations - dependend on kx, ky
	double epsilonk;

public:
	double E[4];
	double dEdlnF[4];
	double dEdlmC[4];
	double dEdlmF[4];
	double dEdq[4];

	double q[4];
	double dqdnF[4];
	double dqdmF[4];
	double dqddF[4];

	Func(const parameters &p, const variables &v) {
		set(p, v);
		//refresh();
	}

	void set(const parameters &p, const variables &v) {
		copyTo(this->p, p);
		copyTo(this->v, v);
		refresh();
	}

	void set(const parameters &p) {
		copyTo(this->p, p);
		refresh();
	}

	// Gutzwiller
	double fq(int s);
	double fdqdnF(int s);
	double fdqdmF(int s);
	double fdqddF(int s);

    // Energy
    double fE(double epsilonk, int s, int l);
    double fdEdlnF(double epsilonk, int s, int l);
    double fdEdlmF(double epsilonk, int s, int l);
    double fdEdq(double epsilonk, int s, int l);

	// REFRESH
	void refresh() {
		nF = v.nF;
		dF = v.dF;
		mF = v.mF;
		mu = v.mu;
		lnF = v.lnF;
		lmF = v.lmF;

		// obliczam funkcje Gutzwillara i poszczegolne pochodne
		for(int i=0; i<2; i++) {
			for(int j=0; j<2; j++) {
				q[j+2*i] 	 = fq(1-2*i);
				dqdnF[j+2*i] = fdqdnF(1-2*i);
				dqdmF[j+2*i] = fdqdmF(1-2*i);
				dqddF[j+2*i] = fdqddF(1-2*i);
			}
		}

		// zeby pozniej mozna bylo wypisac
		glob.qup =  q[0];
		glob.qdown = q[2];

	}

	//REFRESH XY
	void refresh(double kx, double ky) {
		epsilonk = 2.*p.t*(cos(kx) + cos(ky)) + 4.*p.tPrim*cos(kx)*cos(ky);

		refresh();

		for(int i=0; i<2; i++) { // spins: up are E[0] E[1]
			for(int j=0; j<2; j++) { // 1,2 family
				E[j+2*i] 	  = fE(epsilonk, 1-2*i, 1-2*j);
				dEdlnF[j+2*i] = fdEdlnF(epsilonk, 1-2*i, 1-2*j);
				dEdlmF[j+2*i] = fdEdlmF(epsilonk, 1-2*i, 1-2*j);
				dEdq[j+2*i]   = fdEdq(epsilonk, 1-2*i, 1-2*j);
			}
		}
	}
};


/////////////////
///////////////// FullSimplify
/////////////////

double Func::fq(int s) {
	return (2.*Power(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))),2))/((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

double Func::fdqdnF(int s) {
	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((0.25*dF)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          (0.25 + 0.75*Power(dF,2) - 0.5*nF - 0.25*mF*s)/
		           Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
		       1.*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		       2.*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double Func::fdqdmF(int s) {
	return ((dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		       Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
		     (4.*(nF + mF*s)*(1 - 0.5*(nF + mF*s))*((-0.25*dF*s)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          (0.25*(1 + Power(dF,2) - nF)*s)/Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) +
		       1.*s*(nF + mF*s)*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s)))) -
		       2.*s*(1 - 0.5*(nF + mF*s))*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
		          Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
		   (Power(nF + mF*s,2)*Power(1 - 0.5*(nF + mF*s),2));
}

double Func::fdqddF(int s) {
	return 0. + (4.*(dF*Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))*
	      (0. - Power(dF,2)/Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        Sqrt(-1.*Power(dF,2) + 0.5*nF - 0.5*mF*s) +
	        (0. - 4*Power(dF,3) + dF*(-2. + 3.*nF + 1.*mF*s))/
	         (2.*Sqrt((1 + Power(dF,2) - nF)*(-Power(dF,2) + 0.5*(nF + mF*s))))))/
	    ((nF + mF*s)*(1 - 0.5*(nF + mF*s)));
}

///////////////////// ENERGY!

double Func::fE(double epsilonk, int s, int l) {
	//double lmC=0; double A=1., B=1.;
	double A=2., B=6./7.;
	//return (epsilonk - lnF - 2*mu + p.epsilonF - (lmC + lmF + 2*p.h)*s -
	//	     l*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/2.;
	return (epsilonk - lnF - 2*mu + p.epsilonF - (lmF + (A + B)*p.h)*s -
		     l*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2)))/2.;
}

double Func::fdEdlnF(double epsilonk, int s, int l) {
	//double lmC=0; double A=1., B=1.;
	double A=2., B=6./7.;
	//return (-(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s)) -
	//	     Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)))/
	//	   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
	return (-(l*(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s)) -
		     Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2)))/
		   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2)));
}

double Func::fdEdlmF(double epsilonk, int s, int l) {
	//double lmC=0; double A=1., B=1.;
	double A=2., B=6./7.;
	//return -(s*(l*(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s) +
	//        Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2))))/
	//   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
	return -(s*(l*(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s) +
	        Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2))))/
	   (2.*Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2)));
}

double Func::fdEdq(double epsilonk, int s, int l) {
	//double lmC=0; double A=1., B=1.;
	double A=2., B=6./7.;

	//return -((l*Power(p.V,2))/Sqrt(4*Power(p.V,2)*fq(s) + Power(epsilonk + lnF - p.epsilonF - lmC*s + lmF*s,2)));
	return -((l*Power(p.V,2))/Sqrt(4*Power(p.V,2)*fq(s) +
		       Power(epsilonk + lnF - p.epsilonF + (lmF - A*p.h + B*p.h)*s,2)));
}


/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	void onePoint(variables &y, int i, int j);

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->set(p);
	}

	void setParVar(parameters &p, variables &v) {
		this->p = p;
		this->v = v;
		f->set(p, v);
	}

	//void run(variables &w);
	void run2(variables &w);
	void run_global(variables &w);

	double give_fq(int s) { return f->fq(s); }
	double give_fdqdnF(int s) { return f->fdqdnF(s); }
	double give_fdqdmF(int s) { return f->fdqdmF(s); }
	double give_fdqddF(int s) { return f->fdqddF(s); }
	double give_fE(double epsilonk, int s, int l) { return f->fE(epsilonk, s, l); }
	double give_fdEdlnF(double epsilonk, int s, int l) { return f->fdEdlnF(epsilonk, s, l); }
	double give_fdEdlmF(double epsilonk, int s, int l) { return f->fdEdlmF(epsilonk, s, l); }
	double give_fdEdq(double epsilonk, int s, int l) { return f->fdEdq(epsilonk, s, l); }

};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	f->refresh(x, y);
	//if(DEEP_DEBUG) { f->printLambda(); }

	double Etan[4] = { (1-tanh(p.beta*f->E[0]/2.))/2., (1-tanh(p.beta*f->E[1]/2.))/2., (1-tanh(p.beta*f->E[2]/2.))/2., (1-tanh(p.beta*f->E[3]/2.))/2.};

	for(int i=0; i<4; i++) {
		w.lnF += Etan[i]*f->dEdq[i]*f->dqdnF[i];
		w.lmF += Etan[i]*f->dEdq[i]*f->dqdmF[i];
		w.nF += Etan[i]*f->dEdlnF[i];
		w.mu += Etan[i];
		w.mF += Etan[i]*f->dEdlmF[i];
		w.dF += Etan[i]*f->dEdq[i]*f->dqddF[i];

		if(i == 0 || i == 1) { w.magnetization += Etan[i]; }
		else { w.magnetization -= Etan[i]; }

        //if(f->E[i]>0) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
        //if(f->E[i]<0) { w.Energy += f->E[i]; }

		if(DEEP_DEBUG) { printf("########### %f, %f; %d: %f | %e, %e\n", x, y, i, Etan[i], f->E[i], p.beta*f->E[i]); }

        // it would be safer have following code - warning: 30 is chosen arbitrarily
		if(p.beta*f->E[i] > 30.) { w.Energy += -exp(-p.beta*f->E[i])/p.beta; } // basically 0
		else if(p.beta*f->E[i] < -30) { w.Energy += f->E[i]; }
		else { w.Energy += -log(1 + exp(-p.beta*f->E[i]))/p.beta; }
	}
	if(DEEP_DEBUG) {
		printf("###--- lnF=%.15e lmF%.15e nF=%.15e mu=%.15e mF=%.15e dF=%.15e || %.15e %.15e %.15e %.15e | %.15e\n",
				w.lnF, w.lmF, w.nF, w.mu, w.mF, w.dF, f->E[0], p.beta*f->E[0], -exp(-p.beta*f->E[0])/p.beta, -log(1 + exp(-p.beta*f->E[0]))/p.beta, Etan[0]);
	}
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //

void Equa::run2(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation
	int max = p.L/2; // attention: L should be even!

	//int iii = 0;

/*
    // SERIAL CASE - OK
    for(int i=-max; i<max; i++) {
		for(int j=-max; j<max; j++) {
			onePoint(temp, i, j);
     	}
    }
    addTo(temp, w, 1);
    clear(temp);
*/
    // SMART WAY - ALSO OK
    // SERIAL CASE
   	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
   	// * edges (4 times)
   	for(int i=1; i<max; i++) {
   		onePoint(temp, i, i);
   		onePoint(temp, max, i);
   		onePoint(temp, i, 0);
   		//iii += 12;
   	}
      addTo(temp, w, 4);
      clear(temp);
      // * interior area (8 times) - SERIAL CASE
      for(int i=1; i<max; i++) {
   		for(int j=1; j<i; j++) {
   			onePoint(temp, i, j); //iii+=8;									//	sleep(1);
        	}
   	}
   	addTo(temp, w, 8);
   	clear(temp);
   	// * missing corners
   	onePoint(w, max, max);
   	onePoint(w, max, 0);
   	onePoint(w, 0, max);
   	onePoint(w, 0, 0);
   	//iii += 4;

    //printf("%d\n", iii);

	// what was out of sumation
	double L2 = (double)(p.L*p.L);
	w.lnF /= L2;
	w.lmF /= L2;
	w.nF /= L2;
	w.mu /= L2;
	w.mF /= L2;
	w.dF /= L2;
    w.Energy /= L2;
    w.magnetization /= L2;

	w.lnF += v.lnF;
	w.lmF += v.lmF;
	w.nF += v.nF;
	w.mu += -p.n;
	w.mF += v.mF;
	w.dF += 2.*p.U*v.dF;

    v.Energy = w.Energy + v.lnF*v.nF + v.lmF*v.mF + p.U*v.dF*v.dF;
    v.magnetization = w.magnetization;

    // DEBUG
    /*
    double kx = 0.3;
    double ky = 0.2;
    double s = 1;
    double l = 1;
    double epsilonk = 2.*p.t*(cos(kx) + cos(ky)) + 4.*p.tPrim*cos(kx)*cos(ky);
    deb.fq = give_fq(s);
    deb.fdqdnF = give_fdqdnF(s);
    deb.fdqdmF = give_fdqdmF(s);
    deb.fdqddF = give_fdqddF(s);
    deb.fE = give_fE(epsilonk, s, l);
    deb.fdEdlnF = give_fdEdlnF(epsilonk, s, l);
    deb.fdEdlmF = give_fdEdlmF(epsilonk, s, l);
    deb.fdEdq = give_fdEdq(epsilonk, s, l);

	//if(w.chiAB != w.chiAB) { w.chiAB = 1000; }
	//if(w.chiS != w.chiS) { w.chiS = 1000; }
	//if(w.chiT != w.chiT) { w.chiT = 1000; }
	//if(w.barDelta != w.barDelta) { w.barDelta = 1000; }
	//if(w.simDelta != w.simDelta) { w.simDelta = 1000; }
	//if(w.mu != w.mu) { w.mu = 1000; }
	//if(w.mAF != w.mAF) { w.mAF = 1000; }
*/

	if(DEBUG2) {
		printf("###+++ lnF=%.15e lmF%.15e nF=%.15e mu=%.15e mF=%.15e dF=%.15e\n",
				w.lnF, w.lmF, w.nF, w.mu, w.mF, w.dF);
	}

    // results
    case_PAM_s1_v2_g67::glob.Energy = v.Energy + v.mu*p.n;
    case_PAM_s1_v2_g67::glob.magnetization = v.magnetization;
	//glob.qt = f->fgt();
}
/*
void Equa::run_global(variables &w)
{
	clear(w);

	// REFRESH

	case_PAM_s1_v2_g67::n = p.n;
	case_PAM_s1_v2_g67::t = p.t;
	case_PAM_s1_v2_g67::tPrim = p.tPrim;
	case_PAM_s1_v2_g67::U = p.U;
	case_PAM_s1_v2_g67::V = p.V;
	case_PAM_s1_v2_g67::h = p.h;
	case_PAM_s1_v2_g67::epsilonF = p.epsilonF;
	case_PAM_s1_v2_g67::beta = p.beta;
	case_PAM_s1_v2_g67::L = p.L;

	case_PAM_s1_v2_g67::lnF = v.lnF;
	case_PAM_s1_v2_g67::lmF = v.lmF;
	case_PAM_s1_v2_g67::nF = v.nF;
	case_PAM_s1_v2_g67::mu = v.mu;
	case_PAM_s1_v2_g67::mF = v.mF;
	case_PAM_s1_v2_g67::dF = v.dF;

	case_PAM_s1_v2_g67::def_n = 0;
	case_PAM_s1_v2_g67::def_m = 0;
	case_PAM_s1_v2_g67::Energy = 0;

	double y0_d=0, y1_nF=0, y2_lnF=0, y3_mu=0, y4_lmF=0, y5_mF=0;

	  for(int i=0; i<L; i++) {
	      for(int j=0; j<L; j++) {
	         for(int s=-1; s<=1; s+=2) {
	            for(int l=-1; l<=1; l+=2) {
	               double kx = (double)i*2.*M_PI/L;
	               double ky = (double)j*2.*M_PI/L;
	               double epsilonk = 2.*t*(cos(kx) + cos(ky)) + 4.*tPrim*cos(kx)*cos(ky);
	               double E = fE(epsilonk, s, l);
	               double fermi = (1-tanh(beta*E/2.))/2.;

	               y0_d +=     fermi*fdEdq(epsilonk, s, l)*fdqddF(s);
	               y1_nF +=    fermi*fdEdlnF(epsilonk, s, l);
	               y2_lnF +=   fermi*fdEdq(epsilonk, s, l)*fdqdnF(s);
	               y3_mu +=    fermi*(-1);
	               y4_lmF +=   fermi*fdEdq(epsilonk, s, l)*fdqdmF(s);
	               y5_mF +=    fermi*fdEdlmF(epsilonk, s, l);

	               def_n += fermi;
	               def_m += fermi*s;

	               if(beta*E > 30.) { Energy += -exp(-beta*E)/beta; } // basically 0
	               if(beta*E < -30) { Energy += E; }
	               //else { Energy += -log(1 + exp(-beta*E))/beta; }

	               //printf("###%4d %4d %3d %3d: %e\n", i, j, s, l, Energy);
	            }
	         }
	      }
	  }
	  //printf("%f ---- ", Energy);

	  y0_d /=     L*L;
	  y1_nF /=    L*L;
	  y2_lnF /=   L*L;
	  y3_mu /=    L*L;
	  y4_lmF /=   L*L;
	  y5_mF /=    L*L;

	  def_n /= L*L;
	  def_m /= L*L;
	  Energy /= L*L;

	  y0_d +=     2.*U*dF;
	  y1_nF +=    nF;
	  y2_lnF +=   lnF;
	  y3_mu +=    n;
	  y4_lmF +=   lmF;
	  y5_mF +=    mF;

	  Energy += lnF*nF + lmF*mF + U*dF*dF + mu*n;

    // results
	w.lnF = y2_lnF;
	w.lmF = y4_lmF;
	w.nF = y1_nF;
	w.mu = y3_mu;
	w.mF = y5_mF;
	w.dF = y0_d;


    case_PAM_s1_v2_g67::glob.Energy = Energy;
    case_PAM_s1_v2_g67::glob.magnetization = def_m;
	//glob.qt = f->fgt();
    //printf("%f\n", Energy);

    return;
}
*/

int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.n = ((Par *) params)->n();
	p.t = ((Par *) params)->t();
	p.tPrim = ((Par *) params)->tPrim();
	p.U = ((Par *) params)->U();
	p.V = ((Par *) params)->V();
	p.h = ((Par *) params)->h();
	p.epsilonF = ((Par *) params)->epsilonF();
	p.beta = ((Par *) params)->beta();
	p.L = ((Par *) params)->L();

	variables v;
	v.lnF = gsl_vector_get (x, 0);
	v.lmF = gsl_vector_get (x, 1);
	v.nF = gsl_vector_get (x, 2);
	v.mu = gsl_vector_get (x, 3);
	v.mF = gsl_vector_get (x, 4);
	v.dF = gsl_vector_get (x, 5);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    //E.run(w);
    //E.run_global(w);
    E.run2(w);

    // It's end of the calculations

	gsl_vector_set (f, 0, w.lnF);
	gsl_vector_set (f, 1, w.lmF);
	gsl_vector_set (f, 2, w.nF);
	gsl_vector_set (f, 3, w.mu);
	gsl_vector_set (f, 4, w.mF);
	gsl_vector_set (f, 5, w.dF);

	if(DEBUG2) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), gsl_vector_get (x, 5));
		printf("###  %24.15e# %25.15e# %24.15e# %24.15e# %24.15e# %24.15e# | E = %24.15e m = %24.15e\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4), gsl_vector_get (f, 5),
			case_PAM_s1_v2_g67::glob.Energy, case_PAM_s1_v2_g67::glob.magnetization);
		printf("\n");
	}

    return GSL_SUCCESS;
}
/*
void ismaximum(Par &inpar, Var &invar) {
	printf("#"); inpar.print_verboten();
	Var initV = invar;

	// Solving
	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var exact(*res.get_V());
	printf("#Wynik: \n#");
	res.print_label();
	root.print();
	printf("\n");
	initV = exact;

	// EXTRA
	//Var exact = invar;
	//initV = exact;

	// Checking
	printf("#is it maximum?\n");
	printf("#          lnF                  lmF                  nF                       mu                      mF                     d\n");

	struct parameters p;
	p.n = inpar.n();
	p.t = inpar.t();
	p.tPrim = inpar.tPrim();
	p.U = inpar.U();
	p.V = inpar.V();
	p.h = inpar.h();
	p.epsilonF = inpar.epsilonF();
	p.beta = inpar.beta();
	p.L = inpar.L();

	variables v; clear(v);
	Equa E(p, v);

	// EXTRA
	//v.lnF = initV.lnF();
	//v.lmF = initV.lmF();
	//v.nF = initV.nF();
	//v.mu = initV.mu();
	//v.mF = initV.mF();
	//v.dF = initV.dF();
	//variables w; //wynik
	//E.setParVar(p,v);
	//E.run2(w);

	double OptE = glob.Energy;
	double OptEnomu = glob.Energy - inpar.n()*initV.mu();
	while(1) {
		fflush(stdout);
		initV = exact; initV.shake(1e-2);
		initV.print();
		v.lnF = initV.lnF();
		v.lmF = initV.lmF();
		v.nF = initV.nF();
		v.mu = initV.mu();
		v.mF = initV.mF();
		v.dF = initV.dF();

		variables w; //wynik
		E.setParVar(p,v);
		E.run2(w);

		double globE = glob.Energy;
		double globEnomu = glob.Energy - p.n*v.mu;

		double dif = globE-OptE;
		double difnomu = globEnomu-OptEnomu;
		printf("%22.15f %22.15f | %13e", glob.Energy, glob.magnetization, dif); if(dif<0) { printf(" M!"); } else { printf("   "); }
		printf(" | (nomu): %13e", difnomu); if(difnomu<0) { printf(" M!"); }
		printf("\n");

		// EXTRA
		if(difnomu < 0) {
			OptE = globE; OptEnomu = globEnomu;
			initV.set_lnF(v.lnF);
			initV.set_lmF(v.lmF);
			initV.set_nF(v.nF);
			initV.set_mu(v.mu);
			initV.set_mF(v.mF);
			initV.set_dF(v.dF);
		}
	}

	printf("#end ismaximum\n\n\n");
	fflush(stdout);
	return;
}
*/
// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.15e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS: n=%f, t=%f, tPrim=%f, U=%f, V=%f, h=%f, epsilonF=%f, beta=%f, L=%d\n", n(), t(), tPrim(), U(), V(), h(), epsilonF(), beta(), L());
}

void Par::print_label(FILE *file) const {
	fprintf(file, "n   t   tPrim   U   V   h   epsilonF  beta   L\n");
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %23.15e", tabD[i]); }
	fprintf(file, "\n");
}

void Var::print_label(FILE *file) const {
	fprintf(file, "            lnF                    lmF                    nF                       mu                      mF                     d");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_PAM_s1_v2_g67::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%22.15f %22.15f %22.15f %22.15f\n", glob.Energy, glob.magnetization, glob.qup, glob.qdown); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "%22.15f %22.15f %22.15f\n", glob.Energy, glob.qup, glob.qdown); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                       E                       m                       qup                     qdown\n");}

void example() {

	printf("\nEXAMPLE PAM_s2\n");

    Par par(N_PAR);
    par.set_n(1.75);
    par.set_t(-1.);
    par.set_tPrim(0.25);
    par.set_U(12);
    par.set_V(-3);
    par.set_h(0.03);
    par.set_epsilonF(-2.);
    par.set_beta(1500.);
    par.set_L(256);
    par.print_verboten();

    Var initV(N_VAR);
    initV.set_lnF(-1.228966727886987e+00);
    initV.set_lmF(-1.005186920182193e-02);
    initV.set_nF(8.769845923439767e-01);
    initV.set_mu(6.919013389476379e-01);
    initV.set_mF(1.414868211846207e-02);
    initV.set_dF(4.508986650050982e-02);

  //  -1.228966727886987e+00 -1.005186920182193e-02 8.769845923439767e-01 6.919013389476379e-01 1.269263515539406e-02 1.414868211846207e-02 4.508986650050982e-02   -3.659455044347832     0.533274638256591     0.528492344170473
  //		       lnF                  lmF                  nF                     mu                     mC                    mF                   d                      E                      qup                    qdown
  //      -2.270226372436098e+00 5.039109800319238e-03 6.590545498294187e-01 4.646748220161435e+00 3.160851326364592e-03 4.280555711252422e-03 1.733426218259425e-01   -6.570304944582494     0.911264785639954     0.910647838886023
  //      PARAMETERS: n=1.353239, t=-1.000000, tPrim=0.307190, U=11.598375, V=-6.433473, h=0.006557, epsilonF=-2.861979, beta=500.000000, L=256

    printf("Starting point: "); initV.print(); printf("\n\n");
    Res res(N_VAR); res.print_label();
    RootGSL root(par, function, res);
    if(root.Solve(initV) == GSL_SUCCESS) { root.print(); }

    printf("Koniec\n");
}

void example_inf() {
    Par initPar_PM(N_PAR);
    initPar_PM.set_n(1.6);
    initPar_PM.set_t(-1.);
    initPar_PM.set_tPrim(0.25);
    initPar_PM.set_U(40.);
    initPar_PM.set_V(-6);
    initPar_PM.set_h(0.05);
    initPar_PM.set_epsilonF(-3.);
    initPar_PM.set_beta(1500.);
    initPar_PM.set_L(512);
    initPar_PM.print_verboten();

    Par initPar_03hm(N_PAR);
    initPar_03hm.set_n(1.6);
    initPar_03hm.set_t(-1.);
    initPar_03hm.set_tPrim(0.25);
    initPar_03hm.set_U(40);
    initPar_03hm.set_V(-3);
    initPar_03hm.set_h(0.03);
    initPar_03hm.set_epsilonF(-3.);
    initPar_03hm.set_beta(1500.);
    initPar_03hm.set_L(512);

	// I'm choosing start point
	Var initV(N_VAR);

	// INFINITE
	Res res(N_VAR);
	RootGSL root(initPar_PM, function, res);
	printf("sdelta      "); res.print_label(); fflush(stdout);
	while(1) {
		initV.set_lnF(18.*frand()-9.);
	    initV.set_lmF(10.*frand()-5.);
	    initV.set_nF(frand());
	    initV.set_mu(18.*frand()-9);
	    initV.set_mF(frand());
	    initV.set_dF(frand());


		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print();  fflush(stdout);

		}
	}
}


// Przygotowanie do wykresu fazowego [h = 0] (1) [n=1.6, eF = -3] [uwaga! dla U=1000, mamy d na poziomie 0.001 ciągle!]
void run_1D_d() {

    // HALF-METAL
    Par initPar_03hm(N_PAR);
    initPar_03hm.set_n(1.6);
    initPar_03hm.set_t(-1.);
    initPar_03hm.set_tPrim(0.25);
    initPar_03hm.set_U(40);
    initPar_03hm.set_V(-3);
    initPar_03hm.set_h(0.03);
    initPar_03hm.set_epsilonF(-3.);
    initPar_03hm.set_beta(1500.);
    initPar_03hm.set_L(512);

    Var initV_03hm(N_VAR);
    initV_03hm.set_lnF(-3.435693147269268e+00);
    initV_03hm.set_lmF(1.195941721600515e+00);
    initV_03hm.set_nF(6.858444633597155e-01);
    initV_03hm.set_mu(-1.461106592098809e+00);
    initV_03hm.set_mF(4.168186499676897e-01);
    initV_03hm.set_dF(6.394533151762449e-02);

    // -3.201867386359178e+00  -2.158571548452816e-01   6.462739179445612e-01  -1.506532153830258e+00  -8.931859724282243e-02   6.672285303224512e-02    -6.465745279820787     -0.076276752728222      0.621954904677769      0.673322379570407

    // PARAMEGNET
    Par initPar_03PM(N_PAR);
    initPar_03PM.set_n(1.6);
    initPar_03PM.set_t(-1.);
    initPar_03PM.set_tPrim(0.25);
    initPar_03PM.set_U(40);
    initPar_03PM.set_V(-3);
    initPar_03PM.set_h(0.03);
    initPar_03PM.set_epsilonF(-3.);
    initPar_03PM.set_beta(1500.);
    initPar_03PM.set_L(512);

    // -3.435693147269268e+00   1.195941721600515e+00   6.858444633597155e-01  -1.461106592098809e+00   4.168186499676897e-01   6.394533151762449e-02    -6.503847877038575      0.400000000000002      0.784124321836348      0.541311831985667

    Var initV_03PM(N_VAR);
    initV_03PM.set_lnF(-3.435693147269268e+00);
    initV_03PM.set_lmF(1.195941721600515e+00);
    initV_03PM.set_nF(6.858444633597155e-01);
    initV_03PM.set_mu(-1.461106592098809e+00);
    initV_03PM.set_mF(4.168186499676897e-01);
    initV_03PM.set_dF(6.394533151762449e-02);

    // I'm choosing start point
    Par initPar(N_PAR); initPar = initPar_03hm;
	Var initV(N_VAR); initV = initV_03hm;

	printf("#dochodze do punktu startowego g67\n");

	// STANDARD
	//automatic_run_n(initPar, initV, 1.6, 1.80, 0.01);
	//automatic_run_V(initPar, initV, -6, -8, 0.1);
	//automatic_run_U(initPar, initV, 40, 20.0, 1);
	automatic_run_h(initPar, initV, 0.03, 0.01, 0.01);
	automatic_run_h(initPar, initV, 0.01, 0.002, 0.002);
	automatic_run_V(initPar, initV, -3, -1.3, 0.01);
	//automatic_run_h(initPar, initV, 0.005, 0.000, 0.0005);
	//automatic_run_tPrim(initPar, initV, 0.25, 0.0, 0.05); printf("\n");
	//automatic_run_V(initPar, initV, -3., -6., 0.01); printf("\n");
	//automatic_run_epsilonF(initPar, initV, -2, -4.0, 0.1); printf("\n");

	printf("#Jestem w punkcie startowym g67\n");
    Par oldPar = initPar;
    Var oldVar = initV;

    // 1D
    //automatic_run_V(initPar, initV, -3, -6., 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -3, 0., 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_V(initPar, initV, -6, 0., 0.01); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.03, 0.2, 0.0005); initPar=oldPar; initV=oldVar;
    //automatic_run_h(initPar, initV, 0.03, 0.000, 0.0005); initPar=oldPar; initV=oldVar;
    automatic_run_h(initPar, initV, 0.002, 0.25, 0.002); initPar=oldPar; initV=oldVar;

    // 03HM
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.001, -3,  0., 0.01); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.3, 0.001, -3,  0., 0.01); initPar=oldPar; initV=oldVar; printf("\n");

    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.001, -3, -6., 0.01); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.3, 0.001, -3, -6., 0.01); initPar=oldPar; initV=oldVar; printf("\n");

    // 03PM
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.001, -8, -5., 0.01); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.2, 0.001, -8, -5., 0.01); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.0, 0.005, -6, -8., 0.1); initPar=oldPar; initV=oldVar; printf("\n");
    //automatic_2Drun_hV(initPar, initV, 0.03, 0.5, 0.005, -6, -8., 0.1); initPar=oldPar; initV=oldVar; printf("\n");

    //automatic_2Drun_hV(initPar, initV, 0.03, 0.2, 0.001, -6, -3., 0.02); initPar=oldPar; initV=oldVar; printf("\n");

	return;
}


/////////////////////////////////////////////////////////////////////////
///*******************************************************************///
/////////////////////////////////////////////////////////////////////////

void automatic_run_n(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_n, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_n(start + ((double)i)*step);
		printf("%f ", inpar.n());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	inpar.set_n(stop);
	invar = initV;

	printf("#end automatic_run_n\n\n\n");
	fflush(stdout);
	return;
}

void automatic_run_tPrim(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_tPrim, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#tPrim       "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_tPrim(start + ((double)i)*step);
		printf("%f ", inpar.tPrim());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_tPrim\n\n\n");
	fflush(stdout);

	inpar.set_tPrim(stop);
	invar = initV;
	return;
}

void automatic_run_U(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_U, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_U(start + ((double)i)*step);
		printf("%f ", inpar.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_U\n\n\n");
	fflush(stdout);

	inpar.set_U(stop);
	invar = initV;
	return;
}

void automatic_run_V(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_V, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_V(start + ((double)i)*step);
		printf("%f ", inpar.V());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_V\n\n\n");
	fflush(stdout);

	inpar.set_V(stop);
	invar = initV;
	return;
}

void automatic_run_h(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_h, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#h           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_h(start + ((double)i)*step);
		printf("%f ", inpar.h());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_h\n\n\n");
	fflush(stdout);

	inpar.set_h(stop);
	invar = initV;
	return;
}

void automatic_run_epsilonF(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_epsilonF, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#epsilonF    "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_epsilonF(start + ((double)i)*step);
		printf("%f ", inpar.epsilonF());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_epsilonF\n\n\n");
	fflush(stdout);

	inpar.set_epsilonF(stop);
	invar = initV;
	return;
}

void automatic_run_beta(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_beta, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#beta        "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_beta(start + ((double)i)*step);
		printf("%f ", inpar.beta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_beta\n\n\n");
	fflush(stdout);

	inpar.set_beta(stop);
	invar = initV;
	return;
}

void automatic_run_L(Par &inpar, Var &invar, double start, double stop, double step) {
	Var initV = invar;

	Res res(N_VAR);
	RootGSL root(inpar, function, res);
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#automatic_run_L, start %f, stop %f, step %f | ", start, stop, step); inpar.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#L           "); res.print_label();
	fflush(stdout);

	initV = old;
	int LOOPS = int((stop-start)/step);
	if(LOOPS < 0) { LOOPS = -LOOPS; step = -step; }
	for(int i=0; i<=LOOPS; i++) {
		inpar.set_L(start + ((double)i)*step);
		printf("%d ", inpar.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		old.set(res.get_V());
		fflush(stdout);
	}
	printf("#end automatic_run_L\n\n\n");
	fflush(stdout);

	inpar.set_L(stop);
	invar = initV;
	return;
}


///2D
void automatic_2Drun_nV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_n(start1);
	par.set_V(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#n      V        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_n(start1 + ((double)i)*step1);
		par.set_V(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			printf("%f %f   ", par.n(), par.V());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_n(stop1); par.set_V(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_Vn(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_V(start1);
	par.set_n(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#V      n        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_V(start1 + ((double)i)*step1);
		par.set_n(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_n(start2 + ((double)j)*step2);
			printf("%f %f   ", par.V(), par.n());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_V(stop1); par.set_n(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

void automatic_2Drun_hV(Par &par, Var &initV, double start1, double stop1, double step1, double start2, double stop2, double step2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);

	par.set_h(start1);
	par.set_V(start2);
	int LOOPS1 = int((stop1-start1)/step1+0.5); if(LOOPS1 < 0) { LOOPS1 = -LOOPS1+1; step1 = -step1; }
	int LOOPS2 = int((stop2-start2)/step2+0.5); if(LOOPS2 < 0) { LOOPS2 = -LOOPS2+1; step2 = -step2; }

	printf("#"); par.print_verboten();
	printf("#init_p   "); initV.print(); printf("\n");
	printf("#h        V        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	for(int i=0; i<=LOOPS1; i++) {
		par.set_h(start1 + ((double)i)*step1);
		par.set_V(start2);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS2; j++) {
			par.set_V(start2 + ((double)j)*step2);
			printf("%f %f   ", par.h(), par.V());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print(); fflush(stdout);
			initV.set(res.get_V());
		}
		printf("\n");
		fflush(stdout);
	}
	par.set_h(stop1); par.set_V(stop2);
	printf("\n#END\n");
	fflush(stdout);

	return;
}

} // end of namespace case2D2t_SGA_tPrim0
